################################################################################
# Created By:Pietryka
# Creation Date:  2016-12-23
# Purpose: fit conditional fixed effects models on CCES panel data
# Questions: mpietryka@fsu.edu
################################################################################


# PREAMBLE =============================================



# LOAD PACKAGES -----------------
library(Amelia)
library(tidyverse)


# LOAD DATA -----------------

# Workspace from 'CCES2-Impute.R'
load("Data/CCES2-Impute.Rdata")



# DEFINE FUNCTIONS  =============================================

library(survival)

# FUNCTION TO FIT MODEL ON IMPUTATED DATA
felogit_imp <- function(the_formula, the_data) {
  lapply(the_data, function(x){
    clogit(formula = formula(paste0(the_formula, "+ strata(caseid)")),
             data = x)
  })
}


# FUNCTION TO COMBINE ESTIMATES USING RUBIN'S RULES
combine_results <- function(results_list) {
  m_imps <- length(results_list)
  var_names <- results_list[[1]]$coefficients  %>% names()
  b_out  <- sapply(results_list, function(x){x$coefficients})  %>%
    matrix(nrow = m_imps, byrow = TRUE, dimnames = list(NULL, var_names))
  se_out <- sapply(results_list, function(x){summary(x)$coefficients[, 3]}) %>%
    matrix(nrow = m_imps, byrow = TRUE, dimnames = list(NULL, var_names))
  fit_out <- map_df(results_list,~broom::glance(.x))
  ests = mi.meld(q = b_out, se = se_out)
  fit = map_df(fit_out, mean)
  list(q.mi = ests$q.mi, se.mi = ests$se.mi, n = fit$n, aic = fit$AIC, bic = fit$BIC)
}

# DEFINE MODELS   =============================================



outcome_names <- c(
  "turnout",
  "attendmeeting",
  "polsign",
  "workcampaign",
  "donatemoney"
)
controls_a <- c("collgrad", "contacted", "factor(famincquart)")
controls_b <- c(controls_a,
  "married",
  "countymove_change",
  "homeowner",
  "student",
  "employfull",
  "churchfreq",
  "newsint"
)



formulae <- data_frame(y = outcome_names)  %>%
  mutate(form1 = paste0(y, " ~ ", "parent + factor(year)"))  %>%
  mutate(form2 = paste0(y, " ~ ", "parent + factor(year) +",
                            paste(controls_a, collapse = " + ")))  %>%
  mutate(form3 = paste0(y, " ~ ", "parent + factor(year) +",
                            paste(controls_b, collapse = " + ")))  %>%
  mutate_at(vars(form1, form2, form3), funs(map(., formula, env = globalenv())))  %>%
  gather(equation, form, -y)






# FIT MODELS  ====================================



mods_alldata     <- list()
combined_alldata <- list()

for (i in seq_len(nrow(formulae))) {
  print(i)
  # fit models on each imputated dataset
  mods_alldata[[i]] <- felogit_imp(the_formula = formulae[i, "form"]$form,
                           the_data = cces_cleaned)
  # combine estimates  using rubin's rules
  combined_alldata[[i]] <- combine_results(mods_alldata[[i]])
}


# COMBINE RESULTS ====================================



cces_results <- formulae  %>%
  bind_cols(data_frame(
    coefs = map(combined_alldata, 1),
    ses   = map(combined_alldata, 2)
  ))  %>%
  mutate(beta = map_dbl(coefs, ~ .x[colnames(.x) == "parent"]))  %>%
  mutate(se = map_dbl(ses, ~ .x[colnames(.x) == "parent"]))  %>%
  mutate(lb   = beta - 1.96 * se)   %>%
  mutate(ub   = beta + 1.96 * se)  %>%
  mutate(lab = recode(equation,
                      form1 = "1. Baseline Model",
                      form2 = "2. Intermediate Model",
                      form3 = "3. Full Model"))  %>%
  mutate(lab = case_when(
    y != "attendmeeting" ~ lab,
    lab == "1. Baseline Model" ~ paste0(lab, ": No Controls"),
    lab == "2. Intermediate Model" ~ lab,
    lab == "3. Full Model" ~ paste0(lab, ": All Controls")))  %>%
  mutate(lab = forcats::fct_rev(lab))  %>%
  mutate(y_lab = recode(y,
                         'turnout' = "Turnout",
                         'attendmeeting' = "Attend Meeting",
                         'workcampaign' = "Work for Campaign",
                         'donatemoney' = "Donate Money",
                         'polsign' = "Display Sign")
         ) %>%
  mutate(dot_type = ifelse( sign(lb) == sign(ub),
                            "statistically significant",
                            "statistically insignificant")
         )



# SAVE ======================================================

save.image("Data/CCES3-Analysis-Fixed-Effects.Rdata")



